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Abstract 

A one-step, a two-step, an abridged, a skeletal and four detailed kinetic 
schemes of hydrogen oxidation have been tested. A new skeletal kinetic 
scheme of hydrogen oxidation has been developed. The CFD calculations 
were carried out using ANSYS CFX software. Ignition delay times and speeds 
of flames were derived from the computational results. The computational 
data obtained using ANSYS CFX and CHEMKIN, and experimental data 
were compared. The precision, reliability, and range of validity of the kinetic 
schemes in CFD simulations were estimated. The impact of kinetic scheme 
on the results of computations was discussed. The relationship between grid 
spacing, timestep, accuracy, and computational cost were analyzed. 
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1 1. Introduction 

2 In the last decade commercial software packages: ANSYS CFX Flu- 

3 ent {2!, Star-CD jsj], etc. are widely used as a tool for solving Computational 

4 Fluid Dynamics (CFD) problems. According to web site top500.org, compu- 

5 tational capabilities are increasing by a factor of 100 every eight years. The 

6 explosive growth of computational power allows to carry out CFD simulations 

7 using more and more complicated physical models on small clusters of com- 

8 puters or even desktop computers (without using high-power expensive com- 

9 puters). The modern CFD simulations can be multicomponent, multiphase 

10 and multidomain. Heat, mass and radiation transfer as well as chemical 

11 processes can be taken into account in calculations. The increased amount 

12 of model assumptions and parameters is an significant source of errors and 

13 faults. 

14 The processes of verification and validation are very important in CFD 

15 ^. They are ground steps in obtaining a numerical solution (Fig. [1]). The 

16 validation should be done prior to the obtaining of the desired numerical res- 

17 ults while the verification should be done prior the validation. Normally, the 

18 whole numerical model, which includes equations of fluid dynamics, equation 

19 of state and the model of turbulence, is already verified by the developer of 

20 the CFD code and the user should verify only its own user defined models. 

21 Our ultimate aim is the modeling of the flow in a rocket combustion chamber. 

22 In this case we cannot rely on the predefined numerical model, but should 

23 use the models which takes into account the specifics of this complicated 

24 problem. Here we are focusing on the usage of the chemical kinetic models of 

25 hydrogen combustion. In most cases the assumption of thin flame (infinitely 
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26 fast chemical reactions) gives reliable results, so there is no actual need to use 

27 the detailed kinetic mechanisms in CFD simulations. However, the assump- 

28 tion of thin flame is not completely satisfied in rocket combustion chamber 

29 where the turbulence is very high. By this reason the model of the chemical 

30 kinetics should be used for the modeling of the combustion in rocket engine, 

31 but before the model should be verified and validated. 

32 In our case the verification can be done through the comparison with 

33 Chemkin [5] which solves a system of kinetic equations. This gives us a chance 

34 to find and eliminate misprints and to prove that numerical parameters, for 

35 example time-step and grid, do not determine the solution. The next step 

36 should be the validation. After entering into a CFD code a chemical kinetic 

37 model became a part of large physical-chemical numerical model. Generally, 

38 kinetic mechanisms are already validated extensively by their authors, but 

39 after the implementation of the chemical kinetic model the CFD numerical 

40 model needs the validation. Of course chemical reactions drive combustion, 

41 but indeed combustion processes depend on heat and mass transfer too. Al- 

42 though turbulence model, equations of state, transport coefficients, chem- 

43 ical kinetic mechanism can be validated separately, the resulting physical- 

44 chemical model needs the final validation as a whole. 

45 Probably the first example of the verification and the validation of hydro- 

46 gen reaction mechanism in CFD simulations is the work by Mani et al. 6|. 

47 A supersonic flow in a constant-area channel was simulated. The employed 

48 kinetic scheme consisted of 8 reactions without the kinetics of the peroxides. 

49 The supersonic combustion of a hydrogen-air mixture at a high temperat- 

50 ure was modeled. The simple kinetic scheme reproduced experimental data 
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51 properly, what is not surprising when the initial temperature is about 1400 

52 K. 

n 

53 In 1988 Jachimowski [7] modeled scramjet combustion using of a hydro- 

54 gen-air reaction mechanism. He carried out the simulations of hydrogen 

55 combustion at parameters related to flight Mach number 8, 16 and 25. The 

56 hydrogen-air reaction mechanism consisted of 33 reactions, where the de- 

57 tailed hydrogen-oxygen kinetics consisted of 20 reactions. In his study it 

58 was shown that chemical kinetics of HO2 is important at all the studied 

59 Mach numbers. Later Eklund and Stouffer carried out the 3-D CFD 

60 simulation of a flow in a supersonic combustor. They tested two kinetic 

61 models: the detailed model by Jachimowski and the model abridged from 

62 the detailed one. This model was obtained by cutting the kinetics of HO2 

63 and H2O2. The new model consists of only 7 reactions within 6 species plus 

64 bath gas N2. This abridged model is very extensively cited and used due 

65 to the low required computational power in the comparison with detailed 

66 kinetic models. 

67 Kumaran and Babu 8| studied the effect of chemical kinetic models on 

68 CFD calculations. They modeled a compressible, turbulent, reacting flow, 

69 which simulates the supersonic combustion of hydrogen in the jet engine of 

70 hypersonic projectile. The idea of their work is to compare the results ob- 

71 tained using the detailed kinetic model with the previous results obtained 

72 using a single step kinetics. In their past study they attributed the differ- 

73 ence between the numerical results and the experimental data to the inad- 

74 equacy of the one-equation turbulence model and the one-step chemistry. 

75 The simulations with the use of the detailed kinetic mechanism predicts 
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higher combustion efficiency than the calculations using the one-step model. 
The comparison of wall static pressure from the experiments and the numer- 
ical simulations showed that the detailed and the single step chemistry give 
similar results and predict well the positions of pressure peaks, but they fail 



1 the numerical simulations were 

3. 



to predict the values of pressure peaks. A 
carried out using software package Fluent 

The use of a detailed mechanism should give a more precise estimation of 
the main thermodynamic parameters: temperature and pressure, and should 
provide the distribution of the intermediates: OH, H2O2, etc. Generally 
detailed kinetic model provides much more information about combustion 
processes, but the use of it costs an additional CPU power. 



In the considered above wor 
anism of Stahl and Warnatz 



i Kumaran and Babu used the kinetic mech- 
9| published in 1985 as the reference detailed 



mechanism. This mechanism became a little bit old after publication of 



works 



10 



ll| in 2002. In these works the refined data on the rate constant 



of reaction R9 



io| 



H + O2 + M — > HO2 + M 



R9 



and the enthalpy of formation of OH ll| were reported. Since radicals OH 



and HO2 play the essential role in hydrogen oxidation, all kinetic mechanisms 
developed before 2002 should be revised. In the current work the several de- 



tailed hyd 



others 



rogen kinetic mechanisms are tested. One is "old" 



13 



15 



12| , while three 



16| were released after 2002. Besides the outdated thermo 



dynamic data and the outdated rate constants "old" mechanism 
reaction: 

H2 + O2 — ^ OH + OH. 



12| has 



RX 
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The usage of this reaction became marginal 



nowadays, for example: it is 



131]- 15| considered below. The "old" 



88 not included into modern mechanisms 

89 mechanism has been tested in order to see the difference from the updated 

90 mechanisms. 

91 Konnov reported recently about "remaining uncertainties in the kinetic 

92 mechanism of hydrogen combustion" in work 13||. He studied the detailed 

93 hydrogen combustion mechanism. Konnov found two groups of the uncer- 

94 tainties. The first group is associated with the set of the chemical reactions 

95 in the hydrogen-oxygen system. Not all of the possible reactions are in- 

96 eluded in the kinetic mechanisms, and the set of reactions varies from one 

97 mechanism to another. Thereby there is no one conventional set of reactions, 

98 which describes combustion of hydrogen comprehensively, and this problem 

99 is still open. The second type of uncertainties relate to the uncertainties in 

100 the rate constants of the employed reactions. Some of them are not well 

101 defined or the experimental data on them are controversial. It should be 

102 noted that developed by Konnov hydrogen kinetic mechanism [l^ has the 

103 slightly different set of reactions from others extensively cited mechanisms 



104 from Princeton University [17l], National University of Ireland, Galway 15| . 
etc. 



105 



Shatalov et al. 



18|] carried out the analysis of several detailed kinetic 



mechanisms of hydrogen combustion: the mechanism by Konnov 



mechanisms by other authors. Shatalov et al. 18| noted that the use of 



and 



reaction RX does not have a sense, because another parallel channel 

H2 + O2 ^ H + HO2 R-10 
106 has the significantly greater rate constant (in 50 and more times). On the 



other hand the use of reaction RX helps a lot to fit experimental data. They 
pointed out also that below 1100 K reaction R9 effects on the ignition delay 
time significantly. 



There are two very recent 
and Stanford universities 16 



lydrogen reaction mechanisms from Princeton 



19| . The both are validated against the latest 



112 experimental data. For our applications (rocket combustion) the experi 



mental data of Burke and coworkers |20| (the same team as [16]) is the 
most interesting. This experimental data represents the measurements of 
the burning velocities in H2/02/He mixtures at pressures 1-25 atm. These 
measurements are the only available data for hydrogen at high pressure. The 
comparisons of these recent mechanisras with the experimental data showed 
that the mechanism by Burke et al. [16] has a better agreement with ex- 



perimental data at high pressure than mechanism 19|]. By this reason only 
mechanism by Burke et al. [16] has been chosen for the tests in this work. 
A study similar to the current work was carried out by Gerlinger et al. 



2l|. The colleagues studied several hydrogen/air reaction mechanisms in- 



cludin g rn ulti-step schemes 

m 



et al. 



22| and one-step mechanism by Marinov 



23| . The study is focused in the application of reaction mechanisms 



in the simulation of supersonic combustion. The mechanisms were validated 
against ignition delay times. In the validation all mechanism showed sim- 
ilar results excluding one-step mechanism {23 1, which missed a non-Arrenius 
behavior of the experimental data. The authors simulated supersonic com- 
bustion with the different mechanism and compared the results with an ex- 
perimental data. They studied the influence of timestep and numerical grid 
as well. The numerical results showed the sensitivity to the timestep and the 
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132 grid density. Finally the authors conclude that one-step mechanism 23|] is 



15| is more precise 



133 not appropriate, while mechanism by O'Conaire et al. 

134 than other multi-step mechanisms. 

135 While the examples of successful verification, validation and application 

136 of hydrogen reaction mechanisms in CFD simulation of supersonic combus- 

137 tion ramjet exist, the problem of the CFD simulation of hydrogen combustion 

138 in rocket engine is not closed. Scramjet is a specific case and the results ob- 

139 tained for supersonic combustion cannot be extended over the case of rocket 

140 engine. Combustion in rocket engine has its own characteristic features: 

141 high pressures (50-250 atm), the wide span of temperatures from 100 K 

142 to 3500 K, the absence of dilutant (nitrogen). In the case of scramjet the 

143 verification and validation can be done by simulating a supersonic combus- 

144 tion directly what is not possible in the case of combustion in rocket engine 

145 yet. In recent work [24] the group of researchers from five research cen- 

146 ters made the CFD simulations of a flow in a combustion chamber. The 

147 each participant of the project modeled the same object using own meth- 

148 odology. It was the sub-scale rocket engine with 1.5 inch inner diameter, 

149 with one co-axial injector. The combustion chamber had an axial symmetry, 

150 which allowed to carry out the comparison of 2D and 3D modeling. The au- 

151 thors compared steady Reynolds- Average Navier-Stockes (RANS), unsteady 

152 Reynolds-Average Navier-Stockes (URANS) and three different Large Eddy 

153 Simulation (LES) models with the experiment. The comparison showed that 

154 all approaches give the noticeably different results and the only in one case 

155 (LES — stochastic reconstruction model) the obtained results were compar- 

156 able with the experimental data. Indeed the most precise modeling results 
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were obtained with the finest mesh of 255 ■ 10^ cells and using the highest 
computational power of 2 million cumulative CPU hours. However at the 
current moment it is not totally clear how the initial model assumptions af- 
fected the accuracy of the final results, and what assumption or parameter 
impaired the other models. Such comparison is very important from the 
practical point of view because the computational cost and precision vary 
strongly from one numerical model to another. 

The performance of the chemical kinetic models of hydrogen oxidation 
in CFD simulations has been estimated in the current work. The aims of 
the work are the ranking of the selected hydrogen kinetic models and the 
development of the verification, validation and ranking procedures. In the 
current work the performance of the kinetic models are assessed using the 
experimental data on hydrogen ignition 25| and hydrogen fiame 26|]. The 



CFD simulations are carried out using complex physical models. In such con- 
ditions the all-round verification is simply obligatory for the CFD modeling, 
while there is no conventional way to verify combustion models as a part of 
the whole physical-chemical model as well as no conventional way to verify 
and validate physical model itself so far. The simulations have a secondary 
aim to estimate the validity region within the space of the computational 
parameters: mesh, computational scheme, time step. The precision (the dif- 
ference between calculations and experiments) and the computational cost 
(required CPU time) were estimated on the each test case. Global reaction 



model 



23| . two-step scheme |3l|, abridged Jachimowski's mode 



180 skeletal mechanism, four detailed hydrogen mechanisms 



12 



13 



15 



a new 



16| have 



been tested. Thus the results of the work should show what chemical kinetic 
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182 scheme should be used, at which parameters the scheme should be used, how 

183 much computing power it is necessary to have for the fulfilment of a task. 

184 The work is the first step before the CFD simulations of the experiments 



185 carried out at our test facility 



29 



30|. 



186 2. Skeletal kinetic model 

187 In our case as well as in other CFD simulations the numerical and the 

188 physical models have a lot of parameters and need the debugging before 

189 getting the final solution. Parameters, which are not connected with the 

190 problem directly, for example timestep, can seriously obstruct the obtaining 

191 of a solution. Detailed kinetic mechanisms make CFD models too heavy for 

192 the debugging. On the other side global reaction models do not include the 

193 minor species and do not allow to do the all-round debugging (validation and 

194 verification of the properties of HO2, H2O2, etc.). 

195 In this work a skeletal kinetic scheme has been developed, which has the 

196 same set of species as detailed hydrogen mechanisms, but the reduced set 

197 of reactions. This light scheme sped up the formulation of the computa- 

198 tional problem. The problem definition requires to perform the certain set of 

199 the calculations. Different meshes and the models of diffusion and thermal 

200 conductivity were tried before getting the final results. The light skeletal 

201 scheme reduced the amount of the expended CPU hours at the preliminary 

202 stage. The new scheme fills the gap between abridged Jachimowski's model 
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2J, w 



lich has 7 reactions and 6 species, and detailed hydrogen mechanisms 
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151 ] (19-21 reactions, 8 species and bath gases) as well as it allows to 



separate the influence of the amount of reactions and species. 
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The^new scheme was developed from the skeletal model by Kreutz and 
Law [32]. Their skeletal kinetic model has 9 unidirectional reactions and 8 
species. Considering II2/O2 system it may assume that the set of 9 species: 
H2, O2, H2O, H, O, OH, IIO2, H2O2 and bath gas is complete and other 
species (O3, OH^, 011*(A), etc.) can play role only in marginal cases. It 
is necessary to note that at high pressures, what is the case of rocket com- 
bustion chambers, the chain branching proceeds via the formation of HO2, 



II2O2 radicals due to the high rates of the recombination processes 33|. For 
example reaction 

H + 02^0H + 0, Rl 

which is the most important in atmospheric hydrogen flames, is suppressed 
by reaction R9 at pressures above 50 bar. The model by Kreutz and Law 32I 
has a 5 times smaller set of reactions as detailed hydrogen mechanisms and 
can adequately predict the ignition delay times and the ignition limits. On 
the other hand the scheme consists of the irreversible reactions, which means 
that the concentrations of species never reach the equilibrium state. The 
afterburning processes are omitted, which is not important during ignition, 
but leads the mispredictions of species profiles. By these reason the reaction 
set was extended by 6 reactions from detailed hydrogen model [3]. The 
reaction of the quadratic recombination of IIO2 radicals 

HO2 + HO2 H2O2 + O2 

206 was substituted by reactions Rll and R13, see Tables [T], [21 Such extension 

207 increases the computational weight of the model, but it increases the ad- 

208 equacy of the model as well. The new added reactions involve the processes 

209 of radical recombinations, which are important in a post flame zone. 

11 



210 3. Calculations 



211 The CFD calculations have been done with the use of the ANSYS CFX 

212 11 solver \}\, which utilizes the Finite Volume Element Method (FVEM). 

213 The meshes have been created using ANSYS ICEM software. The choice of 

214 the software is given by an adherence to the compatibility of the computer 

215 data and the design documentation. 

216 Two types of tests (simulations) have been done in the work. The first 

217 test case is a quasi 0-dimensional simulation of hydrogen ignition to verify 

218 and validate the models against the experimental data on ignition delay times 

219 25| . The second test case is an 1-dimensional simulation of hydrogen flame 

220 propagation to test the models against the data on the speeds of laminar 



flame 



26|. 



222 An ignition in a perfect adiabatic constant volume reactor has been 

223 modeled as the quasi OD problem. By the formulation the problem is di- 

224 mensionless, but by the settings of calculations it is 3D. The computational 
226 domain represents eighth part of the 1 mm sphere with rigid adiabatic walls. 

226 The mesh consists of 21 nodes and 38 tetrahedron elements. At the ini- 

227 tial moment the whole domain is filled with a stoichiometric hydrogen-air 

228 (O.79N2+O.2IO2) mixture at pressure of 1 atm and temperature in the range 

229 of 900-1400 K. The problem has been solved as a transient task, i.e. the time 

230 evolution of gas conditions has been sought. The laminar model (unsteady 

231 Navier-Stokes equations) was employed, but also one series of simulations 

232 were carried out using the k-e turbulence model. The comparison of results 

233 shows that both models (the laminar and the k-e model) give the same results 

234 in this task. The task was imposed in such way that the results should be 
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235 independent on the choice of turbulence model. Indeed the stagnant homo- 

236 geneous gas mixture is surrounded with the adiabatic rigid walls, so the gas 

237 inside the sphere should be stagnant all time. The object of these calcula- 

238 tions is the estimation of the ignition delay times and the comparison of the 

239 calculated delay times with the data from shock tube experiments 25|, see 



240 Fig. [21 In the calculations the ignition delay times were defined as the time 

241 of a temperature increase up to 500 K relative to the initial temperature. 

242 During the ID tests a freely propagating hydrogen flame has been modeled. 

243 The computational domain consists of 1604 nodes and 400 rectangular prism 

244 elements. All elements are placed along one axis so that the thickness of 

245 the domain equals to one element in two other coordinate axes. The mesh 

246 spacing equals to 5 fim in the direction of the flame propagation. The sep- 

247 arate study of the influence of the grid spacing was carried out where the 

248 spacing was varied from 0.2 fim to 200 fim. The domain represents the rect- 

249 angular with symmetry boundary conditions on the side walls. The domain 

250 has one inlet and one outlet (on the side opposite to inlet). At the outlet 

251 static pressure is specified and equals to 1 atm. At the inlet a hydrogen-air 

252 (O.79N2+O.2IO2) mixture at 298 K and 1 atm flows inside the domain. The 

253 velocity of the mixture is specified at the inlet in the range of 0.5-3.5 m/s so 

254 that the velocity of the flame front reaches a small value in the laboratory 

255 system of coordinates. The mixture composition was varied from equival- 

256 ence ratio of ER = 0.5 to ER = 4.5. The simulations were run as a transient 

257 task. A stationary burning velocity was sought. The laminar model was 

258 employed, also one series of tests were carried out using the k-e turbulence 

259 model. Speed of flame depends essentially on the transport properties of gas, 
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260 SO the temperature dependent thermal conductivity and diffusion coefficients 

261 were used. Thermal diffusion was not taken into account. The system of the 

262 governing equations in ANSYS CFX does not assume mass ffuxes caused by 

263 temperature gradients. 



264 4. Results and discussion 

265 Verification. Comparison with CHEMKIN. 



51]. The results of the 



The both tasks were also solved in CHEMKIN II 
simulations with the help of CHEMKIN II were used as a reference data. 
CHEMKIN is very widely used for solving chemical kinetic problems, where 
the computational problem is formulated as solving of a system of ordinary 
differential equations. Indeed ASYS CFX allows to specify the properties 
of a system by the different ways, while in CHEMKIN task is set in the 
one prescribed format. CHEMKIN uses the modified Arrenius form for rate 
coefficients: 

k = A-T"" -expi-Ea/RT). 

266 The thermodynamic functions: enthalpy, entropy and heat capacity are cal- 

267 culated using the NASA polynomial forms in CHEMKIN. During the cal- 

268 culations in ANSYS CFX the same equations were employed for the rate 

269 constants, the thermodynamic functions and the equations of states, so the 

270 comparisons of the results obtained using the different software have a sense. 

271 The results of the comparisons is depicted in Fig. [3] and Fig. HJ where ANSYS 

272 CFX shows the agreement with CHEMKIN. 

273 Using CHEMKIN the ignition delay times were calculated in the assump- 

274 tions of constant volume and adiabatic walls. In this case the problem defin- 

14 



275 itions (the sets of boundary conditions and kinetic equations) correspond 

276 to each other in ANSYS CFX and CHEMKIN. As a consequence the res- 

277 ults of the simulations using ANSYS CFX agree fully with the computa- 

278 tions in CHEMKIN, see Fig. |3l Indeed it is necessary to note that CFX 

279 solves the 3- dimensional Navier-Stokes equations while CHEMKIN uses the 

280 0-dimensional equation of energy conservation. 

281 The next step is the modeling of freely propagating laminar flame. Zeldovich- 

282 Frank-Kamenetskii equation, which connects flame velocity and reactivity, 

283 gives us a clear view on the problem: 




284 where r is the chemical time scale in reaction zone, and a is the coefficient 

285 of temperature conductivity, which summarizes the effect of diffusion and 

286 heat conductivity through Lewis number Le = 1. In contrast to the previ- 

287 ous case kinetic and transport properties have an equal importance in flame 

288 propagation. 

289 The flame speeds were estimated using ANSYS CFX and PREMIX 

290 Fig. m PREMIX is a subroutine of the CHEMKIN which computes species 

291 and temperature proflles in steady-state laminar flames. The transport prop- 

292 erties were estimated using TRANFIT: the another part of the CHEMKIN 

293 collection. Thermal conductivity, viscosity and diffusion coefficients are es- 

294 timated from the parameters of the Lennard- Jones potential and the dipole 

295 moment of species. The same temperature depended coefficients of thermal 

296 conductivity, viscosity and binary diffusion were used in PREMIX and AN- 

297 SYS CFX, but diffusion ffuxes in multicomponent mixture were approxim- 

298 ated by a different way. 
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299 The specific of H2-O2 system is tliat tlie properties of hydrogen (the 

300 hghtest gas) distinguishes strongly from the properties of other species within 

301 the system. By default ANSYS CFX estimates the transport properties in 

302 the mixture of gases by an inappropriate way for H2-O2 system, where the 

303 influence of the fuel on the transport properties is important. It calculates the 

304 coefficient of thermal conductivity and viscosity of gas mixture using a mass 

305 averaging, and the coefficients of diffusion are calculated from the mixture 

306 bulk viscosity. The problem becomes significant in the case of combustion 

307 in rocket engine where the mixture is not diluted by nitrogen. This problem 

308 can be resolved in ANSYS CFX using CFX Expression Language (CEL) and 

309 setting all coefficients by user as it was done in this work. 

310 In H2-O2 system the diffusion coefficients vary from one component to 

311 other in ~6 times: D02/DH ~ a/ {P02 / f^n) ~ "\/32- Even in the simplest 

312 case, where only II2, O2, H2O are taken into account, the gas mixture can 

313 not be assume as a binary mixture or as a solution of light gas in heavy gas 

314 due to the high fractions of H2 and H2O and the differences in the diffusion 

315 coefficients. In the current work the diffusion coefficients are calculated by 

316 the empirical formula 

317 where Wi is the mass fraction of i-species; Xj is the mole fraction of j-species; 



318 Dij is the binary diffusion coefficient 



34| . After that the diffusion coefficients 



319 of individual species are put into the equation which is responsible for the 

320 transport in CFX: 

PiiUmix - Ui) = — —T^-, (3) 

321 where PiiJJmix — Ui) is the relative mass flux of i-species. The equation is 
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322 not solved for an one constraint component (in our case nitrogen), which 

323 mass fraction is calculated from the constraint that the sum of mass frac- 

324 tions of all species equals to 1. PREMIX (CHEMKIN) uses a more accur- 

325 ate definitions of the diffusion and the thermal conductivity in gas mixture, 

326 and takes into account the thermal diffusion of H and H2. There are two 

327 options: "mixture-averaged properties" and "multicomponent properties" 

328 in PREMIX. "Mixture-averaged" option, which was used here, employ eq. 

329 but does not have a constraint species and employs an additional term 

330 — correction velocity, which makes the net species diffusion flux equal to 

331 zero. "Multicomponent" formulation uses the method described by Dixon- 

332 Lewis 281] , where the coefficients are computed from the solution of a system 



333 of equations defined by the L matrix. 

334 The flame velocities obtained with the use of CFX and CHEMKIN coin- 

335 cide practically with each other. The difference in the results, which is small 

336 (Fig. Hj), should be associated with the distinction between the formulations 

337 of the diffusion fluxes. Coffee and Heimerl compared various methods of 

338 approximating the transport properties of premixed laminar flames, in par- 

339 ticular the methods which have been used in CFX and CHEMKIN. They 

340 found that the difference in flame speed is small for these methods, but 

341 the method, which is employed in CHEMKIN, is more accurate than the 

342 method with constrained species (CFX), which is inaccurate in computing 

343 the diffusion velocity for constrained species. As for the comparison with 

344 experimental data it was shown in recent work [36] that such small over- 

345 shooting around the stoichiometry, which is observed in Fig. |H results from 

346 the neglecting Soret effect (thermodiffusion). 
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347 In our case (laminar hydrogen-air flame at 1 atm) the propagation of the 

348 flame is supported essentially by the diffusion H2 into flame zone and the 

349 diffusion of H into preflame zone. That is why it is important to estimate 

350 accurately the diffusion term, in Fig. [5] we can see the effect of the transport 

351 properties on laminar flame. The maximum of the flame speed is shifted to 

352 the higher equivalent ratios where the diffusive fluxes of H and H2 are higher. 

353 4-^- Validation and testing 

354 Let us consider the results of the first "OD" test case, which is depicted in 



355 Fig. [21 The detailed models 



12 



356 while non-detailed kinetic models 



14-16 



agree with experimental data well, 



3l|, abridged Jachimowski's model 



357 [22] and the new skeletal model have the agreement with experimental data 

358 only in the limited range. The kinetic model by Konnov 13| agrees with 

359 experimental data better than other models. In Fig. [2] it is possible to see 

360 the transition from high-temperature kinetics to low-temperature kinetics 

361 around 950 K. Generally the models show the common trend: more details 

362 — higher accuracy. This conclusion is supported by the results of the ID 

363 test case too. It is possible to conclude from the results that one or two 

364 reactions are not enough to describe the ignition of hydrogen. Probably the 

365 sophistication of abridged Jachimowski's model (7 reactions and H, O, OH as 

366 intermediates) is a reasonable minimum for the modeling of hydrogen com- 

367 bustion in the high temperature region (T > 1000 K). For the modeling in a 

368 wide temperature range the formation of H2O2 and HO2 should be taken into 

369 account. Of course it would be very surprising to see that reduced or global 

370 mechanisms can describe ignition in a wide range of parameters. They are 

371 deduced from detailed mechanisms by neglecting the marginal processes, so 
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372 they can not describe the behavior near margins. In most cases the oxidation 

373 of hydrogen proceeds via the formation of HO2 radical, which is not included 

374 in the reduced mechanisms. In the new mechanism the kinetics of IIO2 is 

375 not comprehensive too. 

376 It is necessary to note that shock tube is not a reactor with adiabatic 

377 solid walls. Due to the boundary layer effect the temperature behind re- 

378 fleeted shock wave slowly increases with time. The discrepancy between the 

379 ideal assumptions and reality arise at large residence times, for the most of 

380 the shock tubes after 1 ms. In Fig. [2l where the ignition delay times are cal- 

381 culated using the boundary conditions of adiabatic solid walls, the detailed 

382 models have a "wrong" trend at low temperatures (large residence times). 

383 At large residence time it is necessary to take into account the real conditions 

384 behind reflected shock wave as it was done in 19|]. In this case the actual 
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16| will better than 



385 agreement between experiments and detailed models 

386 on the graph. 

387 The performed simulations give more information about the evolution of 

388 the system than simply ignition delay times. The "classical" behavior was 

389 observed without anything unusual in the all "OD" tests (by this reason it is 

390 not included in the article). All the kinetic models predict similar temper- 

391 ature (or pressure) time-resolved profiles, which have the induction period, 

392 the following temperature (or pressure) rise, which ends with asymptotic be- 

393 havior. The gas temperature (or pressure) of the combustion products is 

394 predicted correctly by the all kinetic models. 

395 In the ID test case the agreement of the simulating data with the ex- 

396 perimental data is better in sum than in the "ignition" case, see Fig. |5l 
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397 Practically the all models agree with the experimental data. The other dis- 



398 tinctive feature of the o 

399 Jachimowski's model 



stained results is the bad agreement of abridged 



23|. 



22j and the good agreement of one-step model 

400 The results of ID simulations can be interpret in terms of eq. ([1]). The sys- 

401 tem has practically the same physical properties in the all ID simulations. 

402 These allow us to conclude that 



U2 V ^1 

403 where indexes 1 and 2 designate the attribute to different kinetic models. 

404 The ignition delay times should be taken at flame temperature. In our case 

405 the flame temperature amounts ~2000 K. Abridged Jachimowski's model 

406 [22^ has the lowest effective activation energy among the models (see Fig. [2]) 

407 and predicts the significantly larger r at high temperatures. As for one-step 



408 model [23|, which predicts the shortest r at high temperatures, it does not 

409 include atomic hydrogen. It means that the assumptions, which lead us to 

410 eq. (jl]), are not correct for this model. In terms of eq. ([T]) the neglecting of 

411 the diffusion of hydrogen atom should lead to the smaller value of a and to 

412 the essentially less flame velocity, but it is compensated in this model by the 

413 small ignition delay time. 

414 In Figures [2] andb |5] we can see the difference between the results of the 

415 simulations and the experiments. While in the case of the global or reduced 

416 kinetic models the discrepancy can be attributed to the weakness of models, 

417 detailed reaction mechanisms 14, 15| represent the state-of-the-art view 



on hydrogen kinetics. The agreement of these models with experimental data 
and the role of intermediates (H, O, OH, IIO2 and II2O2) were discussed in 
details in orig 
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421 via the formation of highly active intermediates. The experimental study 

422 of the kinetics of the intermediates of H2/O2 system has some difficulties 

423 at temperatures below 900 K where in experiments it is necessary to keep 

424 constant conditions during long residence times. Thus even a detailed kinetic 

425 scheme can fail near ignition and fiammability limits. 

426 Numerical parameters such as time step, grid spacing, type of difference 

427 scheme, etc. should not determine the results of modeling. The proper values 

428 of time step and grid spacing should correspond via the coefficients of physical 

429 model to physical time and space scales. On practical ground the upper limits 

430 of time step and grid spacing are more important, because computational cost 

431 is generally inversely proportional to timestep and the amount of nodes (for 

432 the employed grid the amount of nodes is reciprocally proportional to the 

433 grid spacing), see Fig. [61 The employed values of time step and grid spacing 

434 are normally close to the upper limit. Generally chemical processes can not 

435 be faster than several fractions of microsecond, and hydrodynamic processes 

436 can not take place on a scale smaller than a micrometer. Thus these values 

437 can be set as a reasonable lower limit for the time step and the grid spacing. 

438 The upper limit is quite specific to the details of a task. It is necessary to 

439 estimate the maximum time and mesh steps in each case separately. At a 

440 too high time step the solution diverges. Numerical noise and residuals can 

441 be used as the measure of the proximity to the upper limit of time step. In 

442 this work the adaptive time step has been used and it has been defined by 

443 a residual. The time step was decreased or increased until the value of the 

444 residuals reached the desired level. To estimate the upper limit of mesh step 

445 several simulations were carried out with different spacing, see Fig. [61 On 
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446 the plot we can see a plateau for the cell size below 5 /zm. The upper limit, 

447 which is located near 10 fim, is related to the flame thickness. In flame front 

448 the concentration of hydrogen increases in 2-4 times each 10 /im. Probably 

449 the maximum grid spacing is universal for all hydrogen-air flames at 1 atm 

450 and is defined by the transport properties of the system while the maximum 

451 time step is individual for each kinetic model. As we can see later, maximum 

452 time step is related to the stiffness of kinetic scheme. 

453 Grid spacing, the physical dimensions of computational domain and com- 

454 putational cost are connected with each other. At these conditions grid spa- 

455 cing can limit the applicability of kinetic model. In Table |3]the computational 

456 costs, which are required for the simulation of the evolution of the system 

457 during 1 ms on 1 CPU (Pentium 4) at 2 GHz, are presented for the all tested 

458 models. In the identical the conditions computational cost varies by orders 

459 from one model to another. It is impossible from the data of Table [3] to see 

460 any direct connection of the computational cost with the number of reac- 

461 tions and species. However there is a trend: detailed kinetic models require 

462 a much more computational power than reduced models. Thus the high com- 

463 putational cost limits the application of detailed kinetic mechanisms in CFD 

464 calculations seriously. The grid spacing has the limit near 40 fim after which 

465 the simulations give the absolutely unrealistic results. Detailed models can 

466 be employed only in a special tasks with the computational domains of small 

467 sizes due to the high computational cost and the small grid spacing. 

The computational cost increases strongly from the global reaction mod- 
els to the detailed kinetic mechanisms, but the number of chemical equations 
and species does not completely determine computation cost. Another es- 
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sential parameter is the stiffness of kinetic model. Stiffness is the embedded 
parameter of each kinetic model. It determines the maximum timestep during 
calculations. It arises due to differences in timescales for different chemical 
reactions. To solve kinetic equations it necessary to integrate equations using 
timestep related the smallest timescale over the time interval related to the 
largest timescale. Stiffness can be characterized the ratio of timescales. For 
example, reaction R9 is in 10^° times faster than reaction 

H + O2 ^ O + OH Rl 

in preflame zone while in postflame zone it is vice versa, reaction Rl is faster 
than reaction R9 in 100 times. It is typical for combustion problems, when 

important fast reactions run on the background of slow equally important 
processes. Another example of the embedded stiffness is reaction 

H2O2 + M — ^ OH + OH + M. R15 

468 This reaction determines concentration of H2O2 in flame, and it changes its 

469 direction from reverse to forward after passing the flame front. 

470 5. Conclusions 

471 The eight different kinetic models of hydrogen oxidation were verified, 

472 validated and tested in the CFD simulations what was done using ANSYS 

473 CFX 11 software. The two cases: ignition in adiabatic constant volume 

474 reactor and the propagation of free laminar flame were considered. The 

475 veriflcation of the kinetic models was done through the comparison with the 

476 results obtained with the help of the Chemkin software. The veriflcation 
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477 allowed to eliminate misentries and to define correctly the thermodynamic, 

478 kinetic and transport properties. 

479 The subsequent validation showed that the detailed kinetic schemes are 

480 more precise than the reduced. While it was not found any dependence 

481 between the "speed" of kinetic model and the number of reaction and species, 

482 the reduced kinetic scheme are faster than detailed. The simulations showed 

483 the common trend for kinetic models: more details — higher computational 

484 cost — higher precision. The simulation of the ignition of hydrogen-air 

485 mixture showed that the results are sensitive to the choice of kinetic model. 

486 However in the case of the flame propagation the results are more sensitive 

487 to the model of the transport properties while the reasonable results can be 

488 achieved even with the use of global reaction mechanism. 



The comparison of the simulating data with the experimental data 251.126 1 



showed that detailed kinetic schemes 



14 



16| agree with experiments well. 



491 while the non-detailed schemes agree with the experiments only within a 

492 limited range. The kinetic model by Konnov [l^ has the best agreement 

493 with the experimental data among the tested models. The application of 

494 reduced kinetic schemes of hydrogen combustion, which do not take into 

495 account chemical reactions with HO2 and H2O2, is possible only with strong 

496 limitations. 

497 For the debugging purposes the new skeletal kinetic scheme was developed 

498 which represents the good compromise between computational cost and ac- 

499 curacy. 

500 The carried out study showed that computational results are affected 

501 by the parameters of physical and numerical models. A large amount of 
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502 model parameters is the potential source of errors. The number of different 

503 coefficients reaches thousands in the simulation with the use of a detailed 

504 kinetic model. The parameters of the model can be verified and validated 

505 using the proposed method. 

506 The application of kinetic models in CFD calculations requires the con- 

507 siderable amount of computational power. The maximum time is limited by 

508 the stiffness of model and alters from model to model while the maximum 

509 grid spacing is more or less universal and defined by the thickness of flame 

510 front. 
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Table 1: New skeletal mechanism. 



No 


Ref. No 


Reaction 


A 


n 




Ref. 


1 


Rl 


H + O2 — ^ OH + 


1. 910+14 


0.0 


16.44 


[32] 


2 


R2 


H2 + ^ H + OH2 


5.08e+4 


2.67 


6.292 


[32] 


3 


R3 


H2 + OH — ^H + HaO 


2.16e+8 


1.51 


3.43 


[32] 


4 


R5 


H2 + M i — > H + H + M 


4.57e+19 


-1.4 


105.1 


[15] 


5 


R6 


+ + M^ — >02 + M 


6.17C+15 


-0.5 


0.0 


[15] 


6 


R7 


H + + M — > OH + M 


4.72e+18 


-1.0 


0.0 


[15] 


7 


R8 


H + OH + M — > H2O + M 


4.5e+22 


-2.0 


0.0 


[15] 


8 


R9 


H + O2 + M — > HO2 + M 


6.17e+9 


-1.42 


0.0 


[32] 


9 


RIO 


H + HO2 ^ H2 + O2 


1.66e+13 


0.0 


0.82 


[15] 


10 


R-10 


H2 + O2 ^ H + HO2 


3.68e+13 


0.203 


54.46 


[32] 


11 


Rll 


H + HO2 — > OH + OH 


1.69e+14 


0.0 


0.87 


[32] 


12 


R13 


OH + HO2 ^H20 + 02 


2.89e+13 


0.0 


-0.5 


[15] 


13 


R15 


H2O2 + M — ^ OH + OH + M 


1.2e+17 


0.0 


45.5 


[32] 


14 


R-17 


H2 + HO2 ^H + H202 


3.42e+12 


0.202 


27.12 


[32] 



k = A - T" • e'iqy{—Ea/RT); units: mol, cm'^, K, kcal; thermodynamic data [5|; the reverse 
rate constants (R5, R6) are calculated from the forward rate constants through the 
equilibrium constants. 



Table 2: Efficiency factors for third body term. 



Ref. No 


H 


H2 


H2O 


H2O2 


H02 





O2 


OH 


R5 


1.0 


2.5 


12 


1.0 


1.0 


1.0 


1.0 


1.0 


R6 


0.83 


2.5 


12 


1.0 


1.0 


0.83 


1.0 


1.0 


R7 


0.75 


2.5 


12 


1.0 


1.0 


0.75 


1.0 


1.0 


R8 


1.0 


0.73 


12 


1.0 


1.0 


1.0 


1.0 


1.0 


R9 


1.0 


2.5 


12 


1.0 


1.0 


1.0 


1.0 


1.0 


R15 


1.0 


2.5 


12 


1.0 


1.0 


1.0 


1.0 


1.0 
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Table 3: Parameters of kinetic models and the computational costs. Case 1 — test case 
"ignition", 38 cells; case 2 — test case "flame propagation", 400 cells. 



Model 


Number 
of 

equations* 


Number of 
species 


Case 1, 

L^r^U IlOUrS/I 


Case 2, 
i£IPTJ hours/r 


Marinov et al. 




f 


3 (H2, O2, 
H2O) + N2 


0.18 


0.14 


Lee and Kim c 


ii 


2 


4 (3 + H) + 
bath gas 


0.43 


0.20 




ibridged 
Jachimowski's 
22 1 


7 


6 (3 + H, 0, 
OH) + bath 
gas 


1.7 


5.2 


Zhukov (this w 


Drk) 


f3 


8 + bath gas 


1.6 


19 


Gutheil et al. 


12| 


21 


8 + bath gas 


2.5 


67 




D'Connaire et al. 

Q 


23 


8 + bath gas 


1.8 


71 


Konnov 14\ 




29 


8 + bath gases 


2.3 


36 


Burke et al. [li 


1 


22 


8 + bath gases 


1.8 


35 



The number of equations could exceed the number of reactions because of the possible 
presence of double reactions and of third-body reactions where the activation energy 
depends on the colhsional partner. 
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Figure 1: Logic scheme of validation and verification. 
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Figure 2: Ignition dela y ti mes of a stoichiometric hydrogen-air mixture at 1 atm. Squares 
— experimental data [25|; 1) dash cyan line — Marinov et al. 23|; 2) short dash dot 



magenta line — Lee and Kim [3]J; 3) dash dot purple line — the abridged Jachimowski's 
model 23; 4) short dash line green — Zhukov (this work); 5) solid red line — O'Conaire 
et al. [l5|; 6) dash dot dot line blue — Gutheil et al. [12]; 7) solid black line (the closest 
to exp. data) — Konnov [l^; 8) short dot orange line — Burke et al. 
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Figure 3: Comparison of the simulating data obtained using ANSYS CFX and CHEMKIN. 
Squares — the experimental ignition delay times of a stoichiometric hydrogen-air mixture 
at 1 atm 25|; the kinetic model by O'Conaire et al. dash line (B-spline) and crosses 
— ANSYS CFX, big circles — CHEMKIN. 
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Figure 4: Comparison of the simulating data obtained using ANSYS CFX and CHEMKIN. 
Squares — the experimental burning velocities of a hydrogen-air mixture at 1 atm and 
298 K 26| ; the kinetic model by Konnov ^U] : dash line (B-spline) and crosses — ANSYS 
CFX, big circles — CHEMKIN. 
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Figure 5: Burning velocities of a hydrogen-air mixture at 1 atm. Squares experimental 
data [26|: 1) dash cyan line — Marinov et al. [23]; 2) short dash dot magenta line — Lee 
3) dash dot purple line — the abridged Jachimowski's model 4) short 



and Kim 
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dash line green — Zhukov (this work); 5) solid red line — O'Conaire et al. [15]; 6) dash 
dot dot line blue — Gutheil et al. [12]; 7) solid black line (the closest to exp. data) — 
Konnov [3]; 8) short dot orange line — Burke et al. [igI ]. 
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Figure 6: Simulating results and computational cost as the function of grid spacing, 
hydrogen-air mixture at 1 atm, kinetic model by Zhukov (this work). 
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